// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

// This is a C++ translation/implementation of the work presented in:
// Attilio Milanes, "Tracking magnetic equipotential curves for general combinations of multipolar fields", 2022, CERN, EDMS Nr: 2792136
// also see the ptmf.py script also published on CERN EDMS system

#ifndef DM_DF_MULTIPOLE_HH
#define DM_DF_MULTIPOLE_HH

#include <armadillo> 
#include <cassert>
#include <memory>

#include "distfun.hh"
#include "dfpolygon.hh"

// code specific to Rat
namespace rat{namespace dm{

	// shared pointer definition
	typedef std::shared_ptr<class DFHarmonic> ShDFHarmonicPr;

	// harmonic
	class DFHarmonic: public cmn::Node{
		// properties
		private:
			arma::uword num_poles_;
			fltp reference_radius_; // [m]
			fltp flux_density_; // [T]
			bool is_skew_;

		// methods
		public: 
			// constructor
			DFHarmonic();
			DFHarmonic(
				const fltp reference_radius,
				const arma::uword num_poles, 
				const fltp flux_density, 
				const bool is_skew = false);

			// factory
			static ShDFHarmonicPr create();
			static ShDFHarmonicPr create(
				const fltp reference_radius,
				const arma::uword num_poles, 
				const fltp flux_density, 
				const bool is_skew = false);

			// setters
			void set_num_poles(const arma::uword num_poles);
			void set_flux_density(const fltp flux_density);
			void set_is_skew(const bool is_skew = true);
			void set_reference_radius(const fltp reference_radius);

			// getters
			arma::uword get_num_poles()const;
			fltp get_flux_density()const;
			bool get_is_skew()const;
			fltp get_reference_radius()const;

			// calculate coefficient
			std::complex<fltp> calc_coefficient()const;

			// validity check
			virtual bool is_valid(const bool enable_throws)const override;

			// serialization
			static std::string get_type(); 
			void serialize(
				Json::Value &js, 
				cmn::SList &list) const override;
			void deserialize(
				const Json::Value &js, cmn::DSList &list, 
				const cmn::NodeFactoryMap &factory_list, 
				const boost::filesystem::path &pth) override;
	};


	// shared pointer definition
	typedef std::shared_ptr<class DFMultipole> ShDFMultipolePr;

	// circle distance function for distmesh
	class DFMultipole: public DistFun{
		// properties
		private:
			// poletype
			bool is_tangent_ = false;

			// tolerance for runge kutta
			fltp tolerance_ = 1e-7;

			// pole width
			fltp pole_width_ = 0.06;

			// harmonic content
			std::map<arma::uword, ShDFHarmonicPr> harmonics_;
			
			// aperture settings
			fltp aperture_a_ = 0.05;
			fltp aperture_b_ = 0.03;
			fltp aperture_n_ = 2;

			// working plane x
			fltp x1_ = -RAT_CONST(0.1);
			fltp x2_ = RAT_CONST(0.1);
			arma::uword num_x_ = 41llu;
				
			// working plane y
			fltp y1_ = -RAT_CONST(0.1);
			fltp y2_ = RAT_CONST(0.1);
			arma::uword num_y_ = 41llu;

		// methods
		public:
			// constructcor
			DFMultipole();
			DFMultipole(const std::list<ShDFHarmonicPr>& harmonics);

			// factory
			static ShDFMultipolePr create();
			static ShDFMultipolePr create(const std::list<ShDFHarmonicPr>& harmonics);


			// editting harmonics
			void add_harmonic(
				const ShDFHarmonicPr& harmonic);

			// function to check if harmonics are dipole only
			bool is_pure_dipole()const;

			// marching squares algorithm
			arma::Col<fltp> vertex_interp(
				const fltp level,
				const arma::Col<fltp> &p1,
				const arma::Col<fltp> &p2,
				const fltp valp1,
				const fltp valp2,
				const fltp tol);
			arma::Mat<fltp> marching_squares(
				const arma::Mat<fltp>&R, 
				const arma::Row<fltp> &val, 
				const fltp level);

			// complex functions
			arma::Mat<std::complex<fltp> > calc_F(
				const arma::Mat<std::complex<fltp> > &z)const;
			std::complex<fltp> calc_F(
				const std::complex<fltp> &z)const; // overloaded scalar version
			arma::Mat<std::complex<fltp> > calc_dFdz(
				const arma::Mat<std::complex<fltp> > &z)const;
			arma::Mat<std::complex<fltp> > calc_d2Fdz2(
				const arma::Mat<std::complex<fltp> > &z)const;

			// some handy cases
			arma::Mat<std::complex<fltp> > calc_dAdz(
				const arma::Mat<std::complex<fltp> >&z)const;
			arma::Mat<std::complex<fltp> > calc_dVdz(
				const arma::Mat<std::complex<fltp> >&z)const;
			arma::Mat<std::complex<fltp> > calc_B(
				const arma::Mat<std::complex<fltp> >&z)const;
			arma::Mat<fltp> calc_Bmag(
				const arma::Mat<std::complex<fltp> >&z)const;

			// // calculate curvature
			// arma::Mat<std::complex<fltp> > calc_n_curvature(
			// 	const arma::Mat<std::complex<fltp> >&z)const;

			// create aperture coordinates
			arma::Row<std::complex<fltp> > create_aperture(
				const arma::Row<fltp> &t)const;
			std::complex<fltp> create_aperture(
				const fltp t)const; // overloaded scalar version

			// create polygons
			ShDFPolygonPr create_poles()const;

			// get bounding
			virtual arma::Mat<fltp>::fixed<2,2> get_bounding() const override;

			// distance function
			virtual arma::Col<fltp> calc_distance(const arma::Mat<fltp> &p) const override;

			// validity check
			virtual bool is_valid(const bool enable_throws)const override;

			// serialization
			static std::string get_type(); 
			void serialize(
				Json::Value &js, 
				cmn::SList &list) const override;
			void deserialize(
				const Json::Value &js, cmn::DSList &list, 
				const cmn::NodeFactoryMap &factory_list, 
				const boost::filesystem::path &pth) override;
	};

}}

#endif
